Martin et al. BMC Systems Biology 201 2, 6:54 
http://www.biomedcentral.eom/1752-0509/6/54 



Systems Biology 



METHODOLOGY ARTICLE Open Access 



Assessment of network perturbation amplitudes 
by applying high-throughput data to causal 
biological networks 

Florian Martin 11 , Ty M Thomson 31 , Alain Sewer 1 " 1 , David A Drubin 3 , Carole Mathis 1 , Dirk Weisensee 2 , Dexter Pratt 3 , 
Julia Hoeng 1 and Manuel C Peitsch 1 



Abstract 

Background: High-throughput measurement technologies produce data sets that have the potential to elucidate 
the biological impact of disease, drug treatment, and environmental agents on humans. The scientific community 
faces an ongoing challenge in the analysis of these rich data sources to more accurately characterize biological 
processes that have been perturbed at the mechanistic level. Here, a new approach is built on previous 
methodologies in which high-throughput data was interpreted using prior biological knowledge of cause and 
effect relationships. These relationships are structured into network models that describe specific biological 
processes, such as inflammatory signaling or cell cycle progression. This enables quantitative assessment of network 
perturbation in response to a given stimulus. 

Results: Four complementary methods were devised to quantify treatment-induced activity changes in processes 
described by network models. In addition, companion statistics were developed to qualify significance and 
specificity of the results. This approach is called Network Perturbation Amplitude (NPA) scoring because the 
amplitudes of treatment-induced perturbations are computed for biological network models. The NPA methods 
were tested on two transcriptomic data sets: normal human bronchial epithelial (NHBE) cells treated with the 
pro-inflammatory signaling mediator TNFa, and HCT1 16 colon cancer cells treated with the CDK cell cycle inhibitor 
R547. Each data set was scored against network models representing different aspects of inflammatory signaling 
and cell cycle progression, and these scores were compared with independent measures of pathway activity in 
NHBE cells to verify the approach. The NPA scoring method successfully quantified the amplitude of TNFa-induced 
perturbation for each network model when compared against NF-kB nuclear localization and cell number. In 
addition, the degree and specificity to which CDK-inhibition affected cell cycle and inflammatory signaling were 
meaningfully determined. 

Conclusions: The NPA scoring method leverages high-throughput measurements and a priori literature-derived 
knowledge in the form of network models to characterize the activity change for a broad collection of biological 
processes at high-resolution. Applications of this framework include comparative assessment of the biological 
impact caused by environmental factors, toxic substances, or drug treatments. 



* Correspondence: Alain.Sewer@pmi.com 
+ Equal contributors 

1 Philip Morris International R&D, Philip Morris Products S.A., Quai Jeanrenaud 
5, Neuchatel, 2000, Switzerland 

Full list of author information is available at the end of the article 

O© 2012 Martin et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BlOlVlGCl Central Commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly cited. 



Martin et al. BMC Systems Biology 201 2, 6:54 
http://www.biomedcentral.eom/1752-0509/6/54 



Page 2 of 18 



Background 

Acquisition of large-scale data sets representing a variety 
of data modalities has become a crucial aspect of experi- 
mental system characterization. This strategy enables the 
broad capture of biological information in a short time 
and with a relative small investment of effort, in the hope 
that valuable biological insights might be gained. How- 
ever, the amount of information collected can be over- 
whelming, making interpretation of the data difficult and 
subsequent detailed biological understanding elusive. Re- 
ducing the complexity of such data by evaluating it in a 
relevant biological context is required in order to gain 
meaningful insight. 

High-throughput measurements can be evaluated against 
literature-curated "cause and effect" relationships extracted 
from the Selventa Knowledgebase (see Methods). As illu- 
strated in Figure la, a structure called a "HYP" (derived 
from "hypothesis") is used. A HYP is a specific type of 
network model comprised of a set of causal relation- 
ships connecting a particular biological activity (e.g., 
the increase in abundance or activation of a particular 
kinase, or a more complex network model describing 
a growth factor signaling pathway) to measurable 
downstream entities (e.g., gene expression values) that it 
positively or negatively regulates. Reverse- causal reasoning 
(RCR) uses the measurable downstream entities of a HYP 
to deduce information about the activity of the upstream 
entity of the HYP, based on an appropriate set of mea- 
surements (e.g., differential gene expression from a 
treatment versus control experiment) [1]. Using mea- 
sured downstream effects to deduce the activity of up- 
stream entities is advantageous in that, for gene 
expression data, it does not depend on the "forward" as- 
sumption that mRNA expression changes are always dir- 
ectly correlated with protein activity changes [2-4], an 
assumption that does not take into account the effects of 
translational or post-translational regulation on protein 
activity. 

The HYP can describe causal relationships between an 
upstream biological activity and any type of high- 
throughput data. However, the work described here fo- 
cuses on the evaluation of whole-genome mRNA expres- 
sion changes; thus, the HYP is the equivalent of a causal 
"gene expression signature" for a given entity or process, 
for example, the activity of a particular kinase. Previous 
work has explored the importance of uncovering a char- 
acteristic signature of gene expression changes that 
results from one or more perturbations to a biological 
process, and the subsequent scoring of that signatures 
presence in additional data sets as a measure of that pro- 
cess's specific activity amplitude. Most work in this field 
has involved identifying and scoring signatures that are 
correlated with a disease phenotype [5-8]. These pheno- 
type-derived signatures provide significant classification 



power, but often lack a mechanistic or causal relationship 
between a single specific perturbation and the signature. 
Consequently, these signatures may represent multiple 
distinct unknown perturbations that lead to, or result 
from, the same disease phenotype. 

Alternatively, a number of studies have focused on 
measuring causal signatures based on very specific up- 
stream perturbations, either performed directly in the 
system of interest [9] or coming from closely-related 
published data [10]. Based on the simple, yet powerful 
premise that modulation of cellular pathways and the 
components therein is associated with distinct signatures 
in downstream measurable entities, causally-derived sig- 
natures hypothesize that the "cause" of the signature can 
be identified with high specificity from the measured "ef- 
fect" [11]. These studies have demonstrated the great po- 
tential of applying a causal-pathway activity scoring 
strategy to clinical problems. For example, they have pro- 
vided prognosis predictions in gastric cancer patients 
and indications of specific drug efficacy [10]. 

As a consequence, coupling specific causal HYPs cap- 
tured in the Selventa Knowledgebase with a measure of 
perturbed activity would be a means to further realize 
this clinical potential, as well as the potential to increase 
basic biological understanding that is harbored within 
high-throughput data. However, the HYP infrastructure 
has been previously exclusively employed as an explora- 
tory tool for identifying relevant perturbed biology by 
drawing qualitative mechanistic inferences based on stat- 
istical enrichments [12-14]. Therefore, new methods are 
required to confer an explicit, more quantifiable estimate 
of the degree of HYP activity for a more quantitative 
comparative assessment infrastructure. Such methods 
scoring network perturbation amplitudes (NPA) would 
facilitate a high-resolution comparison of biological 
states, both by virtue of a continuous scale of scores and 
the breadth of HYPs that are immediately available for 
scoring. 

To assess HYP amplitude, four complementary scoring 
algorithms were developed: Strength, Geometric Perturb- 
ation Index (GPI), Measured Abundance Signal Score 
(MASS), and Expected Perturbation Index (EPI). NPA 
scoring was then applied to different inflammatory and 
cell cycle-related HYPs using two transcriptomic data 
sets: a TNFa dose and time series in normal human 
bronchial epithelial (NHBE) cells and a CDK inhibitor 
R547 dose and time series in HCT116 colon cancer cells 
[15]. This study establishes the use of a broad, literature- 
derived knowledgebase to score the amplitude of various 
aspects of biology, which can be defined as very specific 
mechanisms (such as an individual protein activity) that 
are directly proximal to the data, or as a larger network 
of interest that is composed of a collection of individual 
mechanisms. 
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Figure 1 Constructing and scoring causal network models, (a) The Selventa Knowledgebase is composed of individual causal relationships 
that are curated from available sources (generally published scientific literature). For example, in the statement "A— >exp(X)" entity A may 
represent the activity of a particular kinase that has been shown to lead to increased expression of gene X. A HYP is a set of causal relationships 
derived from the Selventa Knowledgebase that relate an upstream entity to the downstream measurable entities that it regulates. In the HYP 
example entity A regulates the expressions of genes W, X, Y, and Z following the specific regulation signs "— >" (positive regulation) and "-|" 
(negative regulation), (b) An aggregated HYP can be generated from a causal network model, which is a collection of causally linked entities 
derived from the Selventa Knowledgebase (grey edges). This causal network model can be augmented with the causal relationships to the 
downstream measureable entities of its nodes (black edges), also derived from the Selventa Knowledgebase. As a consequence, a causal network 
model can be equivalently viewed as a collection of causally linked HYPs. Next, all the downstream measurable entities of the network model 
nodes are combined by adjusting their signs based on the causal relationships in the network. This essential step toward the construction of the 
aggregated HYP is well-defined as long as the network model is causally consistent (see Methods). Because these nodes have a negative causal 
relationship with the reference node (node A), the regulation signs "— >" or "-|" of the downstream measureable entities from nodes C and D are 
inverted (red edges) when constructing the aggregated HYP from the causal network model, (c) NPA scores are calculated from high-throughput 
data, such as differential gene expression data obtained from treatment versus control comparisons, and applied to an NPA scoring algorithm in 
the context of a specific causal network model represented by its aggregated HYP. 



Martin et al. BMC Systems Biology 201 2, 6:54 
http://www.biomedcentral.eom/1752-0509/6/54 



Page 4 of 18 



Results 

The HYP is the foundation for scoring network models 

The HYP represents the relationships between a set of 
measured data, here gene expression data, and a bio- 
logical entity that is a known controller of those genes. 
Additionally, these relationships include the sign (posi- 
tive or negative) of influence between the upstream 
entity and the differential expression of the downstream 
genes. The downstream genes of a HYP are drawn from 
a database of literature-curated causal biological know- 
ledge (Figure la). The causal relationships of a HYP that 
link the upstream entity to downstream genes are the 
substrate for the calculation of process amplitude using 
the NPA scoring algorithms (see below). 

A more general causal network model can be con- 
structed from a set of HYPs that are themselves causally 
connected by literature-derived edges (see Methods). 
Such a network model can be thought of as providing 
higher-level connections between HYPs by linking the 
upstream controllers of these HYPs, based on the path- 
ways graph structure. Complex biological processes such 
as cell proliferation or cellular stress can be efficiently 
described by causal network models [16,17]. 

A complex causal network model of biological entities 
can be transformed into a single HYP by collecting the 
individual HYPs representing entities in the model and 
regrouping the connections of all the downstream genes 
to a single upstream process representing the whole 
complex causal network model; this in essence is a flat- 
tening of the underlying graph structure (Figure lb). In 
this fashion, the activity changes of the biological entities 
described by the network model can be assessed via the 
aggregation of its individual HYPs, such that the under- 
lying gene expression measurements contribute to the 
network as a whole (see Methods for a detailed descrip- 
tion of how the resulting aggregated HYP is 
constructed). 

Scoring HYPs with four NPA methods 

NPA scoring applies a defined algorithm to an experi- 
mental data set consisting in a series of treatment versus 
control comparisons, where the experimental data is 
filtered down to a particular scope of biology (and thus a 
particular set of gene expression relationships) by the 
context of a defined biological network model 
(Figure lc). Specifically, a series of NPA metrics were 
developed to evaluate the activity of the biological 
entities represented by a HYP. The NPA metrics were 
designed such that positive values mean increased activ- 
ity of the biological entities represented by the HYP 
(compared to control), and negative values mean 
decreased activity (compared to control). Furthermore a 
positive or negative relative difference between two NPA 
scores denotes the same relative difference in the 



magnitude of the activity of the biological entities repre- 
sented by the HYP. 

In this study, gene expression data was used to demon- 
strate the NPA approach using four different scoring 
methods: Strength, GPI, MASS, and EPI (see Methods). 
Briefly, Strength is the mean treatment-induced differen- 
tial expressions of the HYPs downstream genes, adjusted 
for the sign of their causal connection to the upstream 
entity of the HYP. GPI is similar to Strength, except that 
the contribution of each gene is additionally weighted by 
taking into account the statistical significance of its dif- 
ferential expression. MASS is the change in absolute 
downstream quantities in a direction supporting an in- 
crease in the upstream entity (i.e., the sum of the mea- 
surements corrected for their causal connection to the 
upstream entity), divided by an average of a total abso- 
lute quantity of the downstreams. Finally, EPI is a 
"smoothed" version of GPI, obtained by averaging a 
slightly modified form of GPI over all possible values of 
a threshold for statistically significant differential gene 
expressions. Each method has specific yet complemen- 
tary advantages (Table 1), having been tailored for a par- 
ticular measurement technology or having different 
conceptual assumptions applied (see Methods). Further 
complementary aspects of these methods will be dis- 
cussed below (see Discussion). 

Statistical annotation of NPA scores 

Each NPA score, regardless of algorithm, represents an 
abstracted view of a set of biological measurements in 
the context of a particular HYP. As such, dozens, hun- 
dreds, or even thousands of measurements may be 
aggregated into a single score. In order to better 
characterize an NPA score and derive value from its use, 
additional statistics that qualify the score are required. 
Two such statistics, Uncertainty and Specificity, were 
developed. Uncertainty is a confidence interval for a par- 
ticular NPA score, while Specificity tests whether an 
NPA score is specific to the downstream genes repre- 
sented by a particular HYP, and not due to a general 
trend of the data (see Methods). 

NPA scoring of an NF-kB HYP accurately assesses NF-kB 
activity 

In order to evaluate the NPA approach, an NF-kB HYP 
was scored for a well-understood and controlled experi- 
mental system - TNFa-treated NHBE cells. The NPA 
results were then compared with an explicit measure of 
the NF-kB complex activity provided by its nuclear 
translocation. 

Activation of the stress- and immune-response tran- 
scription factor NF-kB (nuclear factor kappa -light-chain 
enhancer of activated B cells) has been well-defined as a 
major mediator of TNFa-induced signaling in a variety 
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Table 1 NPA Method Characteristics 


Method 


Main features 


Assumptions 


Pros 


Cons 


Strength 


Linear, unbiased. 

Based on log 2 differential 
measurements (e.g., log 2 fold 
changes of measurements). 


The contributions from the noisy 
downstream measurables sum 
up to zero. 


Intuitive. 


Noisy/biased signals can 
artificially decrease/increase 
the results. 


GPI 


Based on log 2 differential 
measurements. 

Down-weights weak differential 
measurements using false 
non-discovery rates. 


The noisy downstream measurables 
have low false non-discovery rates 
which can be used to minimize 
their contributions. 


Intuitive. 

More robust to noisy 
signals than Strength. 


False non-discovery rate depends 
on the number of experimental 
replicates. 


MASS 


Linear and unbiased in absolute 
non-log 2 scale. 

Dependent on absolute 
changes in measurements. 


Absolute changes in measurements 
are more important than relative 
changes. 


Intuitive. 


Measurements must be directly 
comparable across all downstream 
measurables. 


EPI 


Based on log 2 differential 
measurements. 

Up-weights strong differential 
measurements without using 
false non-discovery rates. 


The downstream measurables with 
higher differential values should 
have stronger contributions than 
those with lower differential values. 


More robust to noisy 
signals than Strength. 

Highest sensitivity to 
strong differential 
measurements. 


Less intuitive. 

Bootstrapping is needed for 
calculating Uncertainty. 



of systems [18,19]. NHBE cells were treated with four 
different doses of TNFa (0.1, 1, 10 and 100 ng/mL) and 
total RNA was collected for microarray measurement at 
four different times after treatment (30 minutes, 2 hours, 
4 hours and 24 hours) (see Methods). All treatments 
were compared to time-matched mock-treated controls 
to obtain 16 contrasts (4 doses x 4 time points). 

Each amplitude scoring method was investigated using a 
HYP created to be a specific measure of NF-kB activation, 
the NF-i<B-direct HYP (Additional file 1). This HYP is 
composed of 155 genes (curated from 247 distinct refer- 
ences, some genes being supported by more than one refer- 
ence) known to be directly regulated by NF-kB (genes 
whose expression is controlled in an NF-i<B-dependent 
manner and whose promoter sequences are directly bound 
by NF-kB). Each scoring method showed the same pattern 
of response to TNFa, having demonstrated a dose- 
dependent response at all times, and a time-dependent 
response that generally saturated at later times (2 hours or 
4 hours, depending on the scoring method; Figure 2a). 
However, there were some differences between the score 
patterns for each scoring method. The closely-related 
Strength and GPI methods produced almost indistinguish- 
able patterns of response, suggesting the contributions 
from noise were balanced in this experiment (see Table 1). 
The EPI method was qualitatively different from Strength 
and GPI in that EPI scores continued to increase from 2 
hours to 4 hours to 24 hours, while Strength and GPI 
scores plateaued from 4 hours to 24 hours. Also, the EPI 
method produced near-zero scores for 0.1 ng/mL TNFa. In 
general, EPI scores for other HYPs and data sets also 
appeared to reduce to 0 (or near to 0) scores that trended 
relatively lower by other methods (Figures 2 and 3). The 



MASS method qualitatively differed from Strength and 
GPI primarily at the 4 and 24 hour time points, with 
Strength and GPI scores increasing from 2 hours to 4 
hours and plateauing from 4 hours to 24 hours, while 
MASS scores plateaued from 2 hours to 4 hours and 
increased from 4 hours to 24 hours. Strength and GPI 
scores met the Specificity criterion (Specificity j^-value < 
0.05) for all conditions. The lowest dose and earliest time 
point for MASS, and the lowest dose for all but the 2 hour 
time point for the EPI method, were found to not be spe- 
cific to the NF-i<B-direct HYP. 

Next, NF-i<B-direct HYP scores were compared to 
NF-kB nuclear translocation. Upon activation, NF-kB is 
transported into the nucleus where it acts to regulate the 
expression of many genes [18,19]. A series of feedback 
loops then lead to the subsequent translocation of NF- 
kB back to the cytoplasm, and this oscillatory cycle con- 
tinues several times [19]. Because NF-kB oscillations 
occur with slightly different periods in different cells in 
the population, the first oscillation is the most reliable 
population-measure of NF-kB activation. Although the 
time of the first oscillation depends on dose, 30 minutes 
after TNFa treatment is a realistic time to measure 
NF-kB nuclear translocation for the doses used here 
(Additional file 2) [19]. Each scoring method produced a 
monotonic, and in some cases nearly linear, relationship 
between score and nuclear translocation, with Pearson 
correlation coefficients between 0.85 and 0.98 for the dif- 
ferent scoring methods (Figure 4). Interestingly, this 
dose-dependent relationship was preserved at different 
times after TNFa treatment (Additional file 3). These 
findings demonstrate that the HYP-based NPA scores 
can quantify NF-kB transcriptional activity. 
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Figure 2 HYP scores for TNFa-treated NHBE cells. Transcriptomic data from TNFa-treated NHBE cells was scored using each scoring method 
(Strength, GPI, EPI, and MASS) for (a) the NF-KB-direct HYP, (b) the IKK/NF-kB signaling HYP, (c) the TNF HYP, and (d) the E2F1 -direct HYP. Error 
bars represent the 95 % confidence interval as determined by the Uncertainty statistic, and scores that failed the Specificity criterion (Specificity p- 
value > 0.05; 1000 comparable HYPs) are shaded in gray. 



NPA scoring of additional HYPs can quantitatively assess 
response to TNFa 

The effects of HYP size and composition were investigated. 
First, the effect of hand-selecting a set of measurements 
that are known to be modulated by NF-kB specifically in a 
TNFa-dependent manner was assessed. A HYP was con- 
structed from a set of 20 genes that were previously mea- 
sured via RT-PCR to assess NF-kB activity in response to 
TNFa treatment in 3T3 mouse fibroblast cells (omitting 2 
genes that have no direct human ortholog) [19]. These 
genes were measured as perturbed by TNFa in 3T3 cells 
upon dosing with TNFa (10 different concentrations span- 
ning 100 ng/mL to 0.005 ng/mL) over a 12 hour time 
course. This HYP produced a very similar pattern of activa- 
tion to the NF-i<B-direct HYP (Additional file 4), suggesting 
that inclusion of genes whose TNFa-dependent expression 
has not been directly verified does not have a detrimental 
effect on the quality of the HYP score. 

Next, the effects of using HYPs derived from upstream 
biological entities that are less proximal to the measurement 
were investigated. To do so, two additional HYPs were con- 
structed: the IKK/NF-kB signaling HYP (Additional file 5), 
which is composed of 992 genes (curated from 414 different 



references) that are known to be modulated by perturbation 
of proteins in a causal network model of signaling from the 
IkB kinase (IKK) proteins to NF-kB activation (Additional 
file 6); and the TNF HYP (Additional file 7), which is com- 
posed of 1741 genes (curated from 589 different references) 
that are known to be modulated by treatment of cells with 
TNFa. The IKK/NF-kB signaling HYP was generated by 
first constructing a causal network model of IKK/NF-kB 
signaling, and then transforming it into a single HYP by the 
aggregation procedure illustrated in Figure lb (see Meth- 
ods). Whereas the NF-i<B-direct HYP is composed entirely 
of genes whose expressions were directly controlled by a 
single transcription factor (NF-kB), each of these two HYPs 
contains genes whose direct transcriptional controller is not 
necessarily known. The expression of these genes may be 
controlled by transcription factors not involved in con- 
structing the HYP. For example, genes in the IKK/NF-kB 
signaling HYP are known to be modulated by perturbation 
of proteins in the IKK/NF-kB signaling causal network 
model, but some of these genes could be regulated as sec- 
ondary effects caused by altered expression of a smaller 
subset of genes that are directly modulated by NF-kB. Also, 
TNFa is a ligand and therefore does not directly mediate 
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Figure 3 HYP scores for CDK inhibitor-treated HCT116 cells. Transcriptomic data from HCT1 16 colon cancer cells treated with the CDK 
inhibitor R547 was scored using each scoring method (Strength, GPI, EPI, and MASS) for (a) the NF-KB-direct HYP, (b) the IKK/NF-kB signaling HYP, 
(c) the TNF HYP, and (d) the E2F1 -direct HYP. Error bars represent the 95 % confidence interval as determined by the Uncertainty statistic. Scores 
that failed the Specificity criterion (Specificity p-value > 0.05; 1000 comparable HYPs) are shaded in gray. 
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the same population of cells. The Pearson correlation between 
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GPI, 0.85 for EPI, and 0.97 for MASS. 



transcription of any genes. Treatment of cells with TNFa 
results in activation of a myriad of transcription factors, any 
of which may directly or indirectly (for example, through 
autocrine signaling) alter the expression of each gene in the 
TNF HYP. 

The IKK/NF-kB signaling HYP and TNF HYP give 
insight into the behaviors of HYPs at different levels of 
proximity to the measurements. The IKK/NF-kB signaling 
HYP is primarily composed of genes that are regulated 
(either directly or indirectly) by NF-kB (Additional file 6), 
and it produced a pattern of response that is very similar 
to the NF-i<B-direct HYP (Figure 2b). This similar pattern 
of response suggests that there is not a large difference 
between the population-level behavior of genes that are 
known to be directly regulated by a transcription factor 
and the behavior of genes where knowledge of direct regu- 
lation is unknown. The time- and dose-dependent 
response that was seen for the NF-i<B-direct HYP appears 
somewhat less robust in the TNF HYP (Figure 2c), for ex- 
ample at the 30 minute time point, but again the four 
methods produced very similar responses. Thus, although 
the general pattern of response was well-preserved among 
the HYPs, minor but noticeable differences in response 
can be observed in HYPs that are less proximal to the 
measurements. 
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NPA scoring detects limited cross-talk between NF-kB and 
cell cycle signaling 

To assess the ability of HYPs to respond specifically to rele- 
vant TNFa signaling perturbations, a HYP was constructed 
for a key cell-cycle component, the transcription factor 
E2F1 (Additional file 8), with the assumption that E2F1 is a 
less direct effector of TNFa signaling compared to NF-kB. 
The E2F1 -direct HYP is composed of 80 genes (curated 
from 54 different references) known to be directly regulated 
by E2F1 (expression controlled by E2F1 and promoter 
sequence bound by E2F1). The E2F1 -direct HYP showed a 
dose-dependent decrease in score for MASS at the 2, 4, and 
24 hour time points, and for Strength and GPI at the 24 
hour time point (Figure 2d). Consistent with this predicted 
decrease in cell cycle progression, CellTiter-Glo® measures 
of cell number found no appreciable increase in cell number 
after 24 hours of TNFa treatment (Additional file 9). Fur- 
ther verification of these conclusions could be performed, 
for example, by measuring the cell cycle progression of the 
sample populations via flow cytometry. 

In order to provide a comparison of NPA results for biol- 
ogy not directly related to NF-kB signaling, the NPA re- 
sponse of the four HYPs introduced above (NF-i<B-direct, 
IKK/NF-kB signaling, TNF, and E2F1 -direct) were assessed 
in response to inhibition of cell cycle progression via a CDK 
inhibitor. Specifically, a publicly available microarray data 
set was used for treatment of HCT116 colon cancer cells 
with three different concentrations of the CDK inhibitor 
R547 (GSE15395) [15] (Figure 3). All four NPA methods 
demonstrated dose- and time-dependent decreases in the 
E2F1 -direct HYP score at the 4 hour, 6 hour, and 24 hour 
time points. The TNF HYP showed a similar pattern of re- 
sponse as the E2F1 -direct HYP, however few of the scores 
passed the Specificity criterion. This suggests that some of 
the genes in this HYP are cell cycle regulated, but are not 
sufficient in number to pass the Specificity criterion. In con- 
trast, the NF-i<B-direct and IKK/NF-kB signaling HYP 
scores did not display this same dose- and time-dependent 
pattern, indicating that these focused HYPs potentially con- 
tain few cell cycle regulated genes. However, the scores for 
the two NF-kB HYPs showed significant increases at the 6 
hour time point (NF-i<B-direct HYP) and the 24 hour time 
point (NF-i<B-direct and IKK/NF-kB signaling HYPs), sug- 
gesting that NF-kB may perhaps be activated by cell cycle 
arrest (for example, [20]). Furthermore, the pattern of scores 
for the NF-i<B-direct and IKK/NF-kB signaling HYPs were 
significantly different for the CDK-inhibition data set, indi- 
cating that these HYPs will not produce near-identical pat- 
terns of scores under circumstances where NF-kB is not 
activated or is perhaps regulated in a complex manner. 

Discussion 

Previous work has demonstrated the utility of exploring 
the reverse causal interpretation of large scale data sets 



using HYPs as opposed to reasoning downstream from 
the data [12-14]. The ability to deduce the degree of 
activity for a broad spectrum of biological processes, 
afforded by an extensive causal knowledgebase, would 
provide enormous potential for facilitating biological 
characterization and yield an even deeper mining of 
information from large-scale data. This approach offers 
the potential to quantify responses of biological systems 
to anything from toxicity and disease processes to thera- 
peutic benefit. This study successfully demonstrated the 
use of the causal connections provided by an appropriate 
knowledgebase as the basis for quantifying the activity 
degree of specific biology from high-throughput data. 
This quantitative application of HYPs, representing pos- 
sibly complex network models, to experimental data 
measuring treatment-induced perturbations is called 
Network Perturbation Amplitude (NPA) scoring. 

Causal directionality is key for the HYP framework 

For all NPA metrics, the proper scoring of a HYP is 
dependent upon the directionality (signs) of the causal 
influences linking the upstream biological entity repre- 
sented by the HYP to the downstream genes whose 
expression it regulates. The knowledgebase harbors in- 
formation about the specific signs (positive or negative) 
of the regulation exerted by the entity represented by 
the HYP on the expression of each of the downstream 
genes. The logic for incorporating differential gene ex- 
pression measurements into an NPA score based on a 
knowledgebase-specified directional blueprint can be 
made via arguments against two specific alternative 
scoring schemes. The case of scoring an activity without 
taking into account the sign of causal influence in the 
HYP can make sense if the HYP represents a transcrip- 
tion factor that always activates or represses genes. 
However, if there are downstream genes in a HYP that 
are controlled in an opposite manner to the others, the 
error of an activity score based on an assumption of a 
single sign becomes apparent: the score contribution of 
genes that are known to be negatively regulated within a 
HYP might cancel, instead of add to, the score contribu- 
tion of genes that are known to be positively regulated 
within a HYP. An alternative tactic would be to 
incorporate the absolute values of the differential 
expression for each gene. This has the problem of 
always producing positive scores for a HYP as well as 
artificially inflating scores: genes that change in a 
manner opposite of how a HYP is known to control 
a gene would add to a HYP activity score rather than 
detract from it. Standard gene enrichment analyses 
usually ignore regulation signs when scoring pathways 
[21,22], but like some newer gene enrichment meth- 
ods [23], NPA assessment methods integrate both types of 
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causal signed relationships of the biology to the measure- 
ments for their output. 

HYPs can be rapidly constructed from an appropriate 
knowledgebase 

HYPs were constructed from the Selventa Knowledgebase, 
a database of causal biological knowledge that allows rapid 
creation of HYPs for any biological process, entity, or caus- 
ally consistent network model that is adequately connected 
to gene expression changes (see Methods). The TNF and 
E2Fl-direct HYPs were created from this knowledgebase 
without any additional literature or experimental investiga- 
tion. For the NF-i<B-direct HYP, because the content of the 
knowledgebase was biologically too limited in this context, 
additional genes that are directly regulated by NF-kB were 
mined from the literature and added to the knowledgebase 
as causal relationships. The additional knowledge concern- 
ing the direct effects of NF-kB was necessary to ensure a 
broadest representation of NF-kB biology. To construct 
the IKK/NF-kB signaling HYP, a network model was built 
by assembling causal relationships between relevant en- 
tities that were represented in the Selventa Knowledgebase. 
Similarly, HYPs and network models can integrate infor- 
mation from other sources besides the Selventa Knowl- 
edgebase, including literature articles and curated 
databases, as long as this information can be interpreted as 
signed causal relationships. The boundaries of HYPs and 
network models are defined during their construction (see 
Methods). For this study the model boundaries were 
chosen based on the biology known to be associated with 
the TNFa treatment of NHBE cells. In a case where the 
expected biology is unknown, a process of identifying biol- 
ogy is required to determine the most appropriate network 
models and HYPs to score. Such an exploratory perspec- 
tive can be provided by evaluating the resulting HYPs 
using the RCR approach [1], which provides a statistical as- 
sessment of whether the activation of a biological entity is 
consistent with measured data, as previously described 
[12-14]. 

Building the IKK/NF-kB signaling network model 
afforded the ability to aggregate the gene expression 
measurements that underlie all the individual HYPS of a 
specific NF-kB network and provide a single score for 
that network. However, in condensing a complex model 
into a single score, there are caveats to consider. Gene 
expressions that have an ambiguous relationship to the 
network (both causally positively and negatively regu- 
lated) must effectively be removed for scoring purposes. 
These ambiguities affect approximately 6 % of the down- 
stream genes of the aggregated IKK/NF-kB signaling 
HYP, which is similar to the case of single HYPs: the NF- 
i<B-direct and the TNF HYPs contains approximately 4 % 
and 7 % ambiguous downstream genes, respectively. 



Their actual impact on the NPA results is expected to re- 
main limited (see below). 

Additionally, resolution with regard to which individual 
entities of the network model are being perturbed is also 
diluted in the overall network score. Thus, when generat- 
ing an aggregated HYP for a network model, key informa- 
tion about the network is not explicitly available, and it is 
important to keep these features in mind when interpret- 
ing NPA scores (see below for a further discussion of the 
methodological perspectives). However, despite these 
caveats, the IKK/NF-kB signaling HYP produced a near- 
identical pattern of response as the NF-i<B-direct HYP, and 
thus is also correlated with the measured physiological 
endpoint, NF-kB activation. 

This finding that similar NPA results were obtained 
using HYPs featuring different characteristics highlights 
an essential aspect of this work: using the same network 
model for calculating the NPA scores of the various ex- 
perimental conditions to be compared (e.g., for all TNFa 
doses and post-treatment times) provides results that are 
robust against having exhaustively captured the per- 
turbed biology in the network model used for NPA scor- 
ing. This aspect is fundamental, especially given the 
practical impossibility of constructing networks models 
capturing all of the biological processes potentially per- 
turbed in a given experiment. It is also exploited when 
constructing network models describing processes that 
are sufficiently generic, e.g., cell proliferation or cellular 
stress [16,17], so that they can be meaningfully applied 
in a variety of experimental situations. 

The robustness of the results also preserves the validity 
of the NPA approach against the possibility of HYPs and 
network models evolving slightly due to the constantly im- 
proving understanding of the biological processes they de- 
scribe. This property was concretely tested with a simple 
step-back calculation consisting of randomly removing 
edges to the HYPs and comparing the corresponding NPA 
results to the original ones. The results demonstrated a 
remarkable robustness: typically, after removal of 20 % of 
their downstream genes, the four HYPs used in this work 
returned GPI profiles that correlated extremely well with 
their original values shown on Figure 2 (Spearman correla- 
tions of 0.99 ± 0.01 obtained on 1000 samples). Therefore 
reasonable future additions to the Selventa Knowledgebase 
are not expected to significantly impact the NPA results 
presented in the work. As a corollary, these robustness 
considerations also support the choice of discarding the 
downstream genes expressions that have an ambiguous re- 
lationship to the HYP upstream entity. Examination of 
Additional file 1, Additional file 5, Additional file 7, Add- 
itional file 8 showed that the fraction of ambiguous down- 
stream genes never exceeds 10 %, which clearly indicates 
that this effect does not affect the NPA results. 
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NPA scoring methods accurately assess biological 
activation 

Four different algorithms were developed to quantify the 
amplitude of perturbation of a HYP. Each method 
employs a different approach to evaluate the degree of 
perturbation between two experimental measures for a 
given HYP (see Table 1 and Methods). Despite their 
differences, the four methods generally produced similar 
qualitative results, suggesting that each NPA scoring 
method is able to effectively quantify the changes in 
activity of the underlying biological processes. This claim 
is supported by the fact that the NPA scores correlated 
well with NF-kB nuclear translocation, a standard meas- 
ure of NF-kB activity (correlation coefficients between 
0.85 and 0.98). This correlation further validates our 
method of constructing HYPs from a database of prior 
knowledge. 

Future work will confirm the circumstances in which 
each method is expected to be the most appropriate. For 
example, the MASS algorithm was developed to use ab- 
solute measurement technologies such as an absolute 
transcript count offered by quantitative next generation 
sequencing. Given more appropriate measurement tech- 
nology, the MASS method may be more applicable in 
circumstances where small differential expressions in a 
set of highly expressed genes have a dominant effect over 
large differential expressions in a set of weakly expressed 
genes. On the other hand, the GPI algorithm, and to an 
even larger extent the EPI algorithm, down-weight the 
contributions from genes with poor statistical signifi- 
cance, favoring small sets of strongly differentially 
expressed genes rather than large sets of weakly differen- 
tially expressed ones. From this point of view, the 
Strength algorithm is unbiased since it contains no 
weighting factors. However, because an NPA score repre- 
sents a condensed view of the biology underlying a HYP, 
the ability to assess the amplitude of its perturbation 
with complementary NPA methods also highlights which 
conclusions are robust versus which conclusions may be 
specific to a particular NPA method. For example, the 
four NPA methods supported the same time- and 
dose- dependent NF-kB activation in response to 
TNFa (Figure 2a), whereas only Strength and GPI 
suggested NF-kB activation in response to CDK inhib- 
ition (Figure 3a). 

There are some important considerations when using 
NPA scoring methods to evaluate HYPs. Scores are 
meant to be directly compared between different treat- 
ment versus control contrasts when using the identical 
HYP. Scores cannot be quantitatively compared between 
two different HYPs without first verifying that the rela- 
tionship between a change in the activity of the HYPs 
upstream entity and the resulting change in the expres- 
sion of downstream genes is conserved between the two 



HYPs. In general, this relationship is not expected to be 
conserved due to differences in the dynamic range of 
expression of individual genes that compose each HYP. 
Additionally, a HYP with a higher number of down- 
stream gene expressions may be expected to represent 
broader biology than a smaller HYP, and thus in any 
given experiment, a smaller fraction of genes in the lar- 
ger HYP may be perturbed, resulting in a lower score 
than a smaller HYP. However, additional statistical 
power is gained in the Uncertainty and Specificity statis- 
tics with increasing number of downstream gene expres- 
sions in a HYP, such that the weaker scores of larger 
HYPs can be just as significant and meaningful as higher 
scores from smaller HYPs. 

Although scores cannot be directly compared between 
two different HYPs, the pattern of scores across a set of 
contrasts can be qualitatively compared. Likewise, the 
absolute magnitude of the NPA scores should not be 
directly compared between two amplitude scoring meth- 
ods, but the pattern of scores across NPA scoring meth- 
ods can be qualitatively compared, keeping in mind that 
the scoring methods may be assessing different aspects 
of the contrasts. 

The NPA score represents an abstracted measure of 
the data represented in the HYP. The score captures the 
amplitude of the perturbation of a HYP, but does not 
capture which genes in the HYP most strongly contrib- 
ute to the score. For example, of the 20 genes that con- 
tribute most to the IKK/NF-kB signaling HYP score 
upon TNFa treatment (100 ng/mL, 24 hour), only one is 
in common with the 20 genes that contribute most to 
the IKK/NF-kB signaling HYP upon CDK inhibition (0.6 
uM, 24 hour). Given that the NF-i<B-direct HYP consists 
of only 155 genes, this suggests that there is a significant 
difference in the biology represented by the NF-kB- 
direct HYP scores in these two cases. 

Uncertainty and specificity of NPA scores 

Uncertainty estimates the confidence interval of each 
NPA statistic, and therefore also tests the nullity of the 
score accounting for the experimental error. The Specifi- 
city statistic gives a measure of whether the score is 
dependent on the expression of specific genes in the 
HYP, or is instead dependent on a particular property 
(the likelihood of modulation) of the ensemble of gene 
expressions in the HYP. Although this definition of Spe- 
cificity is useful, there are some important points to 
ensure that Specificity is interpreted appropriately. First, 
a weak Specificity does not mean that the score fails to 
accurately characterize the amplitude of the process 
described by the HYP. Rather, it means that many other 
comparable HYPs would obtain a similar score. For 
example, a very weak score (approximately zero) for a 
transcriptomic HYP is likely to have a weak Specificity 
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because the majority of the genes on a microarray are 
unchanged under most conditions. Thus, any random as- 
sortment of genes in a HYP might produce a similarly 
low score. Weak Specificity for low scores could there- 
fore be an indication that the genes in the HYP are not 
sufficiently perturbed. Alternatively, a high score with a 
weak Specificity does not indicate that the process mea- 
sured by the HYP is not perturbed. Rather, it indicates 
that comparable gene expressions are perturbed to a 
similar extent, suggesting that other processes with com- 
parable HYPs are likewise perturbed, and thus the score 
cannot be attributed with high probability to the process 
represented in the HYP. For example, the fact that the 
pattern of Strength scores for the TNF HYP in the CDK 
inhibitor experiment is similar to the pattern of Strength 
scores for the E2F1 -direct HYP suggests that the TNF 
HYP may contain some genes that are cell cycle con- 
trolled (Figure 3). However, this number of genes is not 
sufficient to distinguish this score from the "background" 
of scores for comparable HYPs, as only one of the fifteen 
TNF HYP Strength scores met the Specificity criterion. 
In fact, 32 genes are common to the TNF HYP and the 
E2F1 -direct HYP, which constitutes more than a third of 
the E2Fl-direct HYP, but only one fortieth of the TNF 
HYP. Methods such as Network Component Analysis 
[24,25] could possibly be adapted to resolve overlaps be- 
tween HYPs by assigning shared gene expressions to the 
most statistically likely HYP, potentially increasing the 
precision of each HYP and modulating score Specificity 
appropriately. 

Together, the Uncertainty and Specificity statistics en- 
able the identification of non-specific and non-significant 
scores in HYPs when scored against unrelated perturba- 
tions. These statistics demonstrate that TNFa treatment of 
NHBE cells only has a significant effect on cell cycle pro- 
gression when the dose is above 0.1 ng/mL, and that this 
effect takes two-to-four hours to appear. Also, these statis- 
tics support the conclusion that some NF-kB- regulated 
genes are upregulated at 6 and 24 hours after CDK-inhib- 
ition in HCT116 cells, but likely not at 4 hours or earlier. 

Potential applications beyond the comparative 
assessment of biological impact 

The NPA approach developed in this work aims at quan- 
tifying the treatment-induced perturbations of the bio- 
logical processes described by causal network models. It 
enables the comparative assessment of the biological im- 
pact from high-throughput data in response to given 
stimuli. However the NPA approach could be also used 
in more exploratory perspectives. For instance, by apply- 
ing NPA scoring to each HYP in a causal network model, 
rather than constructing and scoring a single aggregated 
HYP for the model (Figure lb), differences in activation 
across a model could be investigated. Scoring individual 



HYPs within a model instead of the larger aggregated 
model HYP presents a tradeoff of increased granularity 
of information at the expense of statistical robustness, 
due to the smaller sizes of the HYPs being scored. An- 
other possibility would be to use the NPA scores and 
their companion statistics to identify which processes are 
potentially activated in response to a given perturbation, 
and thus help guide the construction of a HYP or of a 
causal network model that capture the relevant per- 
turbed biology. 

Finally, NPA scores could be used as a supplementary 
source of information in studying different types of math- 
ematical models of regulatory networks. The fact that NPA 
scores provide quantitative measurements for the response 
of entities that are not explicitly measured or measureable 
can be exploited in the construction, calibration, or evalu- 
ation phases of such models. For instance, in the case of 
TNFa-treated NHBE cells considered in this study, the 
NPA scores provide direct quantitative measurements of 
the inflammatory response of the system, a quantity that 
would be difficult to access in the absence of the NF-kB 
nuclear translocation measurements performed in this 
work. 

Conclusions 

NPA is an integrated approach that combines high- 
throughput experimental data and a knowledge-driven 
HYP, which provides measurable quantities causally 
affected by a targeted biological process, to quantify the 
activity changes of that process relative to a control 
(non-perturbed) state of the system. The utility of the 
NPA method lies in the synergy of on-demand HYP gen- 
eration from an extensive causal knowledgebase with a 
continuous measure of its activity change. Four NPA 
scoring methods with complementary strengths per- 
formed similarly in this study, but individually have the 
potential to wield distinct advantages for specific 
circumstances. Additionally, qualifying NPA statistics 
enabled effective use of these scoring metrics and can be 
applied to similar methods developed elsewhere. When 
applied to TNFa- and CDK inhibitor- treated cell micro - 
array data, NPA scores for NF-kB and cell cycle networks 
correlated with expected dose response relationships and 
specific measured pathway outputs. NPA scoring also sug- 
gested possible cross-talk between NF-kB activation and 
the cell cycle that could be investigated experimentally. 
With a broad spectrum of biology available to score within 
the Selventa Knowledgebase, NPA metrics and statistics 
can be used to assess amplitude of perturbation on many 
orders, from a single molecule to that of a complex, 
higher-order causal network model representing complex 
biological processes. This approach enables a quantitative, 
systems-wide understanding of the biological mechanisms 
leading to diseases. This is the first step towards the 
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development of computational tools designed to compara- 
tively measure any perturbation, including exposure to 
toxic substances, the effects of drug treatment, or patient 
stratification by individual biology. 

Methods 

Experimental procedure 

Normal human bronchial epithelial cells (Lonza Walkers- 
ville, Inc.) were cultured in standard growth medium (Clo- 
netics medium, Lonza Walkersville, Inc.). Cells were either 
treated with TNFa (Sigma) or a vehicle control (HBSS), 
and then harvested after the desired treatment length 
(30 minutes, 2 hours, 4 hours, or 24 hours). Cells were 
immediately put on ice and split into three technical repli- 
cates from which total RNA was extracted using RNeasy™ 
Microkit (Qiagen). The processed RNA samples are then 
hybridized to Affymetrix U133 Plus 2.0 microarrays. Cell 
viability and cell counts were controlled for all conditions 
after 24 hours with CellTiter-Glo® assay (Pro mega). NF-kB 
nuclear translocation was measured using Cellomics 
NF-kB Activation HCS Reagent Kit (Thermo Scientific). 

Data processing and algorithm implementation 

Data processing and NPA methods were implemented in 
the R statistical environment [26]. Raw RNA expression 
data was analyzed using the affy and limma packages of 
the Bioconductor suite of microarray analysis tools avail- 
able in the R statistical environment [27,28]. Robust 
Microarray Analysis (RMA) background correction and 
quantile normalization were used to generate probe set 
expression values [29] a . An overall linear model was fit 
to the data for all groups of replicates, and specific con- 
trasts of interest (comparisons of "treated" and "control" 
conditions) were evaluated to generate raw ^-values for 
each probe set on the expression array. Raw ^-values 
were subsequently corrected for multiple testing effects 
using Benjamini-Hochberg false discovery rate (FDR), as 
described hereinafter. 

Probe sets were matched to RNA Abundance nodes in 
the Selventa Knowledgebase (see below) using the HG- 
U133_Plus_2.na30 probe set mappings and the following 
criteria. First, only "at" or "s_at" probe sets were consid- 
ered. Second, probe sets that mapped to multiple genes 
were discarded. Third, when multiple probe sets mapped 
to the same gene, preference was given to "at" probe sets 
over "s_at" probe sets. Finally, when there still remained 
multiple probe sets mapped to the same gene, the probe 
set with the lowest geometric mean FDR-corrected 
j?-value across all contrasts of interest was selected. A 
linear model was then re-fit for all groups of replicates 
to only those probe sets that mapped to RNA Abun- 
dance nodes in the knowledgebase, and FDR-corrected 
^-values were recomputed as described in this section. 



The Selventa Knowledgebase 

The Selventa Knowledgebase is a comprehensive reposi- 
tory containing over 1.5 million nodes (biological con- 
cepts and entities) and over 7.5 million edges (assertions 
about causal and non-causal relationships between 
nodes). The assertions in the knowledgebase are derived 
from peer-reviewed scientific literature as well as other 
public and proprietary databases (Figure la). Specifically, 
each assertion describes an individual experimental 
observation from an experiment performed in a human, 
mouse, and rat species context, either in vitro or in vivo. 
Assertions also capture information about the referring 
source (e.g. the PubMed ID (PMID) for journal articles 
listed in MEDLINE), as well as key contextual informa- 
tion including the species (human, mouse, or rat) and 
the tissue or cell line from which the experimental 
observation was derived. An example causal assertion is 
the increased transcriptional activity of NF-kB (nuclear 
factor kappa-light-chain-enhancer of activated B cells) 
causes an increase in the mRNA expression of CXCL1 
(Chemokine (C-X-C motif) ligand 1) [HeLa cell line; 
Human; PMID 16414985]. The knowledgebase contains 
causal relationships derived from healthy tissues and dis- 
ease areas such as inflammation, metabolic diseases, car- 
diovascular injury, liver injury, and cancer. 

Constructing HYPs and causal networks models 

In addition to edges that enable the construction of 
HYPs (i.e., causal relationships between one upstream 
entity and its downstream measurables; Figure la), the 
Selventa Knowledgebase also contains literature-curated 
relationships between the upstream entities of different 
HYPs. These edges can be assembled to construct larger 
causal network models describing more exhaustively the 
biological processes under consideration (Figure lb). 

The principles guiding the construction of the HYPs 
used in this work were pragmatic (e.g., the NF-i<B-direct, 
TNF, and E2Fl-direct HYPs, see Results). The down- 
stream measurable nodes derived from the Selventa 
Knowledgebase in an automated manner are gene 
expressions ("exp (...)" nodes on Figure la) that have 
been shown to be causally linked to the HYP upstream 
entity in a variety of experiments, cell types, and even 
species. Ideally, this context information stored in the 
knowledgebase could be leveraged to construct HYPs 
using only knowledge derived from the cell type and per- 
turbation under current study. However, despite the 
large amount of information in the Selventa Knowledge- 
base, there was generally insufficient material from simi- 
lar contexts to enable automated construction of 
context-specific HYPs. Therefore context-based HYP 
filters were not used in this work. This choice is sup- 
ported by the fact that HYPs with a larger number of 
downstream genes increase the likelihood that a number 
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of the downstream genes in the HYP will be relevant to 
the system being investigated. Furthermore, the positive 
and negative score contributions from irrelevant down- 
stream genes are expected to balance each other and 
produce a background noise averaging to zero and thus 
not to affect the final NPA results. This assumption was 
confirmed a posteriori by the comparisons of the 
Strength and GPI NPA scores (see Results). 

Causal network models (e.g. the IKK/NF-kB signaling 
network model, see Results) are constructed by manually 
assembling causal relationships connecting HYP up- 
stream entities derived from the Selventa Knowledgebase 
(Figure lb). While context information is not taken into 
account during the process of HYP construction, the 
process of network model construction is guided by the 
cell type and perturbation relevance of each causal con- 
nection in order to get a high-quality outcome. In this 
study the boundaries of the IKK/NF-kB signaling net- 
work model were fixed so that signaling in NHBE cells 
in response to TNFa treatment was accurately described. 
A more detailed discussion on the construction of causal 
network models is given under the construction process 
of the "Literature Model" in two published studies 
involving similar network types [16,17]. 

It is important to note the difference between the 
knowledge-based causal network models used in this 
work and the networks constructed from expression data 
using network inference approaches [30]. Although both 
are "causal", their content is fundamentally different. 
Knowledge-based network models are constituted of lit- 
erature-curated causal relationships between biological 
entities, for which the NPA approach enables "backward" 
deduction of the activity changes (or perturbation) from 
differential gene expression data. The expression data- 
based inferred networks describe gene-gene interactions 
that are "forward" deduced from the expression values of 
the corresponding genes and measured over a large col- 
lection of experimental conditions. Therefore, whereas 
linked genes are interacting in the networks inferred 
from expression data, gene expression nodes (i.e. "exp 
(...)" nodes in Figure 1) that are connected via a HYP in 
a knowledge-based network model are under the influ- 
ence of a common upstream entity (e.g. the transcrip- 
tional activity of NF-kB for the NF-i<B-direct HYP). 

Constructing a HYP from a causal network model 

A causal network model is composed of multiple caus- 
ally-linked nodes that are biologically related, including 
HYPs [16,17]. In order to generate the corresponding 
aggregated HYP, a reference node must first be selected 
within the network model. The reference node can be 
any entity in the network whose level or activity is posi- 
tively related to the activity of the network as a whole (as 
opposed to, for example, an inhibitor whose activity may 



be negatively related to the network activity). Next, the 
causal relationship between each node in the model and 
the reference node is determined. This can only be done 
by first requiring that the model be "causally consistent" 
(see below). The signs of regulation of downstream 
measurable entities (here, gene expressions) for each 
node in the model are adjusted based on the relationship 
between that model node and the reference node. For 
example, the signs of the downstream gene expressions 
for a model node that has a positive causal relationship 
with the reference node (i.e., that node is expected to be 
positively regulated when the reference node increases) 
are maintained. On the other hand, the signs of the 
downstream gene expressions for a model node with a 
negative causal relationship with the reference node (i.e., 
that node is expected to be negatively regulated when 
the reference node increases) are inverted. All the down- 
stream gene expressions and their signs are then 
assembled into a single HYP (Figure lb), and down- 
stream gene expressions with contradictory signs (from 
multiple model nodes) are omitted from the aggregated 
HYP. 

For a network model to be causally consistent, for 
an increase in any node in the model, it must be pos- 
sible to unambiguously map a sign of "positive regula- 
tion" or "negative regulation" on every other node in 
the model by following the causal relationships that 
connect the nodes. For example, any model with 
negative feedback loops cannot be used to construct 
an aggregated HYP. Biological interpretation can be 
used to resolve ambiguities to construct causally con- 
sistent models by considering what process is being 
scored by the HYP, and in what sign each node is ef- 
fectively related to the reference node. For example, 
the node where a negative feedback connects back to 
the model has a particular relationship with the 
process being scored, and although the negative feed- 
back may regulate this node, it should not change this 
relationship. Thus, the connection between the nega- 
tive feedback loop and this node can be removed 
from the model to obtain causal consistency in a 
manner that is congruent with our biological expecta- 
tions. Such a biology-driven procedure to build caus- 
ally consistent network models implies a careful 
definition of the network model boundaries. Resolving 
the potential sign ambiguities in its nodes involves 
additional refinements in delimiting the biological 
processes described by the model. 

NPA scoring algorithms 
Strength 

The Strength NPA scoring method was developed as the 
simplest measure of HYP amplitude — the weighted 
mean of the measurement differences (here, gene log 2 
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differential expressions) where the weighting factors are 
the signs of the measurable downstream entities in the 
HYP: 

Strength = —Z? =1 Si -fr 

where is the log 2 differential expression (i.e., the log 2 fold 
change) of the i th gene in the HYP, and s t is the sign (+1 for 
positive regulation and -1 for negative regulation) of the i th 
gene in the HYP, and N is the number of measurable down- 
stream gene expressions in the HYP. This translates to the 
sum of log 2 differential gene expressions with positive regu- 
lation in the HYP minus the sum of the log 2 differential 
gene expressions with negative regulation in the HYP, 
divided by the total number of gene expressions in the 
HYP. Thus, a positive Strength score indicates that the 
HYPs downstream gene expressions and their signs of 
regulation are matched within the data, and the process 
described by the HYP is upregulated in the treated condi- 
tion compared to the control condition. A zero Strength 
score indicates that the process is unchanged, and a 
negative Strength score indicates that the process is 
downregulated. 

The Strength scoring method uses data for all mea- 
sured downstream gene expressions in the HYP, regard- 
less of data quality. However, the differential expression 
values used by the Strength method may be dominated 
by noise for low absolute measures of the control condi- 
tion, and thus may provide unreliable data (which should 
be evidenced by high uncertainty). The Strength method 
assumes that this noise is evenly distributed and that 
there are a sufficient number of measurable down- 
streams in the HYP such that measurement noise is aver- 
aged out across all downstreams. 

Geometric perturbation index 

A HYP can be seen as a unit sign vector s = (1,1,-1,1,. ..,-1)/ 
V/V in the TV-dimensional downstream measurable space 
(where each dimension represents a downstream meas- 
urable, here gene expression, of the HYP). The observed 
effect of perturbation on the downstream gene expressions 
is also a vector in this space. So geometrically, the ampli- 
tude of the perturbation in the HYP can be quantified by 
projecting the differential log 2 expression vector onto 
the hypothesis unit vector s. However, the downstream 
measurements of a HYP come from a generic model. 
To deal explicitly with the specificity of data supporting 
an NPA score, each downstream is assigned a belief of 
activation, which is set to be the local false non-discovery 
rate (fndr t - 1-fdr^) [31]. The false discovery rates fdr t are 
obtained from the raw ^-values using the Benjamini- 
Hochberg multiple testing corrections (see above). It is 
equivalent to weight the dimensions of the downstream 
gene expression space according to the belief of each 



differential expression and therefore consider a weighted 
scalar product to define the geometry of the gene ex- 
pression space: <s\P>fndr - s T -diag(fndr)-p. Hence, the 
Geometric Perturbation Index (GPI) scoring method is 
defined as: 

GPI^-^Zl.srJkdn-Pi 

Note that Strength and GPI are closely related - GPI 
is normalized by VN rather than N. By weighting the dif- 
ferential log 2 expression with false non-discovery rate, 
individual differential expression values for which there 
is little confidence are moved closer to zero (no change), 
while values for which there is stronger confidence are 
minimally decreased. A positive GPI score indicates an 
upregulation of the process described by the HYP, a zero 
GPI score indicates that the process is unchanged along 
the direction s of the HYP, and a negative GPI score indi- 
cates that the process is downregulated. 

Measured abundance signal score 

Both Strength and GPI quantify log 2 differential values 
(i.e., log 2 fold changes) of measurements in the HYP. 
However, there may be cases where an absolute change 
in mRNA, protein, or some other measurable physical 
quantity is a better measure of the biological effect on 
the process represented by a HYP. For example, an in- 
crease from 1 to 10 copies of an mRNA transcript (an 
absolute increase of 9 transcripts) may be less significant 
than an increase from 10 to 100 copies of the same tran- 
script (an absolute increase of 90 transcripts). An NPA 
scoring method was devised to quantify absolute changes 
in entities that represent physical quantities, Measured 
Abundance Signal Score (MASS): 

. . , Z+-iSi ' (treatedi-controli) 

MASS = — ^ — — -f- 

LiLi [treated i + control )) /2 

where treated t is the measurement (not in log 2 scale) for 
the i th downstream measurable (here, gene expression) 
in the treated sample, and control is the measurement 
(not in log 2 scale) for the i th downstream gene expres- 
sion in the control sample. The numerator of MASS 
represents the change in the absolute downstream gene 
expression quantities in a direction supporting an 
increase in the process described by the HYP. The 
denominator of MASS represents the average of the total 
absolute quantity of the downstream gene expressions. 
Thus, MASS can be thought of as quantifying the abso- 
lute change in the downstream gene expressions (cor- 
rected for the predicted sign s t of each downstream in 
the HYP) compared to the total quantity of the down- 
streams. Rather than use the total quantity of the down- 
streams in either the treated or control condition alone, 



Martin et al. BMC Systems Biology 201 2, 6:54 
http://www.biomedcentral.eom/1752-0509/6/54 



Page 15 of 18 



the average of the treated and control conditions was 
used to ensure that the MASS scoring method is sym- 
metric about the experimental contrast (i.e., MASS (treated 
versus control) = -MASS(control versus treated)). 

The MASS method is applicable to any measurement 
technique that quantifies physical measurables in a man- 
ner such that measurements are proportional to absolute 
quantities across all measurable entities (i.e., the mea- 
surements for different entities can be compared dir- 
ectly). It should be noted that for microarrays, this 
implementation is a proof of concept as there is only a 
loose correlation between the signal for different probe 
sets and the absolute amount of the transcripts [29,32]. 

Expected perturbation index 

The Expected Perturbation Index (EPI) can be thought 
of as a smoothed version of the GPI. The rationale is to 
consider a HYP as a random variable over a compact 
interval [-M,M] containing all possible differential ex- 
pression values (in log 2 scale, typically M=15 as gene 
expression signal saturates at 15). The downstream en- 
tities (here, gene expressions) of the HYP are used as evi- 
dence to construct the distribution Phyp of the HYP: for 
any 0 in [-M,M], the density of the HYP at 0 is propor- 
tional to the total evidence of "correct predictions" s?fii 
above 0. More precisely: 



^HYP(C|)) 0C < 



ly JAI 

jsj^-'i such that Si-faxp ^ 

ly N 

^y-^j such that 5j-^<(|) j^j- 



if (p > e 
if (p < -8 



The normalization of the measure PHYp(<p)'d(p is done 
through continuous interpolation of p H yp between [-e,e], 
for e sufficiently small (typically on the order of machine 
precision). The Expected Perturbation Index is then 
defined as being the expectation of the HYP with respect 
to the distribution Phyp defined above: 

1 rM 

The EPI expresses the fact that high correctly predicted 
differential expressions s$ t provide stronger evidence for 
the HYP than low ones and will therefore receive higher 
weights. The implementation is done by using the method 
of rectangle on the differential expression values, leading 
to the following formula to compute the EPI: 

EPI = 14 SU ch that (^) f >o(^)r(^ 

o(*0r(^ 2?i(-^)-((-^-(-^) w ) 



+1* 



such that (s/?).<0V 



HYP, respectively, subscripts i refers to the ascendant 
ordering of the and fi 0 = (9. 

Note that the EPI is, in absolute value, bounded by M, 
and that the false non-discovery rates are not used. Indeed, 
high differential expression value are taken into account 
more often than lower ones, and this additional weighting 
to high differential expression enables a more specific meas- 
ure of activity. 

NPA scoring statistics 
Uncertainty 

Each NPA score is a random variable. As such, the statis- 
tical significance of the NPA scores can be assessed by esti- 
mating confidence intervals around the score, or, 
equivalently, by evaluating the null hypothesis of the score 
equaling zero at a given type-I error risk a (often, and 
herein, a = 0.05). In the case of Strength and GPI, this task 
can be completed analytically, while in the case of MASS 
and EPI, a bootstrapping approach is necessary. 

The theoretical distribution of the differential log 2 expres- 
sion is deduced from the estimation and test procedure 
used. In the case of t-statistics or moderated t-statistics 
(e.g., produced from a linear model by the limma R pack- 
age), the (theoretical) distribution of the /?/s is assumed to 
be normal with variance sdf (having df t degrees of freedom). 

In this context, Strength becomes a random variable 
consisting of a weighted sum of independent (approxi- 
mately) normal distributions. As a consequence, the distri- 
bution of the Strength statistic is (approximately) a normal 
distribution itself, with variance Sd 2 Strength = (1/N 2 ) ♦ E sdf. 
Hence, a t-statistic t = Strength/Sd strength can De derived, 
whose degrees of freedom Df are estimated with the 
Welch-Satterthwaite equation [33,34]. Therefore a (1- 
a)-confidence interval (e.g., 95 % confidence interval) for 
the Strength is given by: 

Strength ± t^ -Sd strength 

In the case of the GPI, the situation is almost identical 
to the Strength, except for the additional dependence on 
the weighting factors fndri. In turn, fndri depends on the 
non-adjusted j^-value, and therefore on the original t-sta- 
tistics for the fa, ti = /3i/sdi. The key step toward a test 
statistic for the GPI is the estimation of the variance of 
GPI. To this end, the variance is estimated through a 
first-order Taylor expansion to obtain: 

where: 



where (sf$)i is s$ b n + and n_ are the number of differential 
log 2 expressions positively and negatively predicted by the 



Yi=fridri + 2 



Pi sdt 
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Using the central limit theorem, an approximate t-statistic 
is derived as t = GPI/Sd GPb whose degrees of freedom Df 
are again estimated by the Welch-Satterthwaite equation. 

In the partial derivatives of the GPI with respect to f} it 
the derivative of fndri = 1- fdr t involves two terms: one 
involving the derivative of the Benjamini-Hochberg ad- 
justment factor (fdr/pi), multiplied by p t (the p- value for 
the differential expression of the i th gene), and the other 
one involving the derivative of p t . The former is assumed 
to be zero, because the adjustment factor (fdr/pi) is im- 
portant only when the ^-values p t are small The latter is 
computed analytically. Hence the (1 -a) -confidence inter- 
val for GPI is given by: 

GPI ± t a J f 2 .Sd GPI 

In the case of MASS and EPI, a parametric bootstrap- 
ping of the distributions of /},• was used, taking advantage 
of the normality of the estimators /?/ deduced from the 
statistical approach used to compute the differential gene 
expressions (see above). As the assumption for the applica- 
tion of percentile bootstrapping seems to be violated for 
MASS and EPI, the confidence interval was estimated 
using the bias-corrected percentile method [35,36]. 

Specificity 

It is important to consider whether the computed NPA 
score is specific to the HYP of interest, or is a general prop- 
erty of the entire data set. For example, a score that indi- 
cated a two-fold increase in a given process holds less 
meaning if all measurements in the entire data set also 
increased two-fold. Thus, the Specificity statistic is com- 
puted as a means to identify scores that can be attributed 
with high probability to the specific biological entity or 
process represented by the HYP. Specificity is computed by 
assessing the likelihood of the following null hypothesis: 
"The amplitude score is not representative of the specific 
HYP, but instead is representative of a general trend in the 
data set that can be measured by any HYP that is compar- 
able to the HYP of interest." The first step to computing 
Specificity is to construct a set of comparable HYPs (see 
below). Next, an amplitude score is computed for each of 
these HYPs using the same data set. Finally, Specificity is 
computed as a two-tailed j?-value by placing the amplitude 
score for the HYP of interest on the distribution of scores 
for the comparable HYPs (Additional file 10). Scores that 
have Specificity ^-values less than 0.05 are considered to be 
scores that can be attributed with high confidence to the 
HYP of interest. 

The key to computing the Specificity statistic is con- 
structing relevant comparable HYPs. A simple comparable 
HYP could be composed of downstream measurable 
entities selected at random from the set of all measurable 
entities to produce a HYP of the same size as the HYP of 



interest. However, for some data sets including gene ex- 
pression microarrays, many of the measurable entities are 
highly unlikely to change under any given circumstance 
(for example, genes whose expression are measured with 
ineffective probe sets). The HYP of interest will likely con- 
tain fewer of these entities than any "comparable" HYP be- 
cause these entities are less likely to be affected by the 
process of interest, and are thus unlikely to be included in 
the HYP of interest. Therefore, a comparable HYP con- 
structed by random sampling from all measurable entities 
will be biased towards showing a weaker (or no) change. 
This makes the amplitude score for the HYP of interest 
appear more unlikely to have occurred by chance given 
the null hypothesis - and thus be more specific - than 
might have otherwise been expected (Additional file 10). 

The large body of causal knowledge available was used 
to construct more relevant comparable HYPs for tran- 
scriptomic data by first identifying the number of up- 
stream controllers (distinct entities upstream of a gene in 
causal statements in the knowledgebase) for each gene in 
the entire data set, including the genes in the HYP. The 
number of upstream controllers reflects the number of dif- 
ferent experiments or perturbations that caused the gene 
to be modulated, and acts as a naive estimate for the likeli- 
hood of each gene being modulated in the current experi- 
ment. For example, a gene that is only regulated under 
very specific circumstances is unlikely to be modulated in 
many data sets that are curated in the knowledgebase. 
Thus, only a few entities that causally regulate the gene 
will exist in the knowledgebase. In contrast, a gene whose 
expression is modulated by a large number of experimen- 
tal perturbations is likely to be modulated in many data 
sets that are represented in the knowledgebase, and thus 
the knowledgebase will likely contain knowledge of many 
entities that causally regulate the gene. Therefore, a com- 
parable HYP is constructed by replacing each gene in the 
original HYP with another measured gene with a similar 
number of upstream controllers. 

In order to accomplish this, all measured genes were 
ranked based on their number of upstream controllers in 
the knowledgebase, and divided into cadres of a fixed 
number of measurables. To avoid having a cadre con- 
taining only a few measurables (for example, when the 
number of measurables is not divisible by the desired 
cadre size), the cadre that contained measurables with 
the fewest number of upstream controllers was allowed 
to have more measurables than the other cadres. For ex- 
ample, a cadre size of 100 measurables was used, so the 
cadre that contained measurables with the fewest num- 
ber of upstream controllers had between 100 and 199 
members. Comparable HYPs were constructed by swap- 
ping each measurable in the original HYP for another 
measurable from within the same cadre. In this manner, 
the comparable HYPs were the same size as the original 
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HYP and contained the same distribution of frequently 
and infrequently modulated measurables. This process 
enabled the construction of alternative HYPs that were 
as comparable as possible to the HYP of interest, given 
the knowledge available in the knowledgebase. 

NPA scores were computed for each comparable HYP, 
and these scores were sorted into ascending order. The 
fraction of scores that were greater than and less than 
the score for the original HYP were counted, and the 
lesser of these two values was doubled to arrive at a two- 
tailed ^-value. This ^-value is called the Specificity. 

Comparable HYPs constructed in this manner provide 
a more stringent test to assess the Specificity of ampli- 
tude scores. Given the benefits offered by this method of 
computing Specificity, it was carried forward for this 
study. 

Endnotes 

a The gene expression data used in this publication have 
been deposited in ArrayExpress and are accessible 
through accession number E-MTAB-1027 [37]. 

Additional files 



model architecture (middle). CHUK, IKBKB, and I KB KG act as inhibitors of 
NFKBIA, NFKBIB, and NFKBIE, which are in turn inhibitors of NFKB1, NFKB2, 
and RELA. The nodes used in the model are listed under each section. 
The nodes in bold represent nodes that have downstream gene 
expression measurables in the knowledgebase, and the number of 
measurables is given in the square brackets (because the same 
downstream may be found under multiple nodes, these 1227 
downstream measurables correspond to 992 unique measurables). The 
notations used in the knowledgebase are as follows: "CHUK P@S" 
represents CHUK phosphorylated at serine (where the residue is given if 
known), "CHUK P@ST" represents CHUK phosphorylated at serine or 
threonine (the exact residue is unknown), "kaof(CHUK)" represents the 
kinase activity of CHUK, "CHUK:IKBKB" represents the complex of CHUK 
and IKBKB proteins, "IkappaB kinase complex Hs" represents an aggregate 
of the various IkB kinases (CHUK, IKBKB, and I KB KG) in Homo sapiens (Hs), 
"degradationof(NFKBIA)" represents the process of NFKBIA degradation, 
and "taof(NFKBI)" represents the transcriptional activity of NFKB1. 

Additional file 7: The TNF HYP. Each row of the table contains a causal 
statement extracted from the Selventa Knowledgebase and describes a 
gene known to be modulated by the TNFa treatment of cells. 

Additional file 8: The E2F1 -direct HYP. Each row of the table contains 
a causal statement extracted from the Selventa Knowledgebase and 
describes the connection between E2F1 and one of its direct target 
genes. 

Additional file 9: CellTiter-Glo® measurements of cell numbers. 

CellTiter-Glo® fluorescence measurements of cell number after 24 hours 
of NHBE cell treatment with various doses of TNFa. Fluorescence intensity 
is reported as a percentage of the mock-treated sample, and error bars 
represent the standard deviation of the fluorescence intensity for three 
different fields of view of the same population of cells. 

Additional file 10: Computing Specificity statistics. The histogram of 
MASS scores for HYPs comparable to the NF-KB-direct HYP (1 ng/mL 
TNFa treatment of NHBE cells for 0.5 hour). The top histogram resulted 
from selecting measurables at random from all measurables (Specificity p- 
value of 0), and the bottom histogram resulted from selecting 
measurables with comparable likelihood of modulation (Specificity 
p-value of 0.072). The solid line indicates the NF-KB-direct HYP MASS 
score. The Specificity p-value was computed by doubling the total 
fraction of counts that were greater than this score (the fraction of counts 
in the red boxes). 
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Additional file 1: The NF-KB-direct HYP. Each row of the table 
contains a causal statement describing the connection between NF-kB 
and one of its direct target genes. 

Additional file 2: TNFa dose-dependent induction of NF-kB nuclear 
translocation, (a) Non-treated NHBE cells show a diffuse cytoplasmic 
staining of NF-kB (left) while after adding TNFa to the culture medium 
(right), the nucleus of stimulated cells is strongly labeled (magnification 
1 0X/0.3); (b) NF-kB nuclear fluorescence intensity (Fl) per dose of TNFa. 
For each group, the nuclear fluorescence intensity was measured in 500 
cells per well of three replicates. Across-group comparisons by one-way 
ANOVA test were all p < 0.001 . 

Additional file 3: HYP scores versus NF-kB nuclear translocation. 

The NF-KB-direct HYP scores for each amplitude scoring method 
(Strength, GPI, EPI and MASS) and each time point (30 minutes, 2 hours, 4 
hours, 24 hours) plotted against NF-kB nuclear translocation at 30 
minutes. Score error bars represent the 95 % confidence interval as 
determined by the Uncertainty statistic. Error bars in NF-kB nuclear 
translocation represent the standard deviation of the mean nuclear 
translocation for three different fields of view of the same population of 
cells. 

Additional file 4: Comparison of NF-KB-direct HYP scores with 20- 
gene NF-kB HYP scores. Transcriptomic data from TNFa-treated NHBE 
cells was scored using each scoring method (Strength, GPI, EPI and 
MASS) for (a) the NF-KB-direct HYP, (b) a HYP composed of 20 NF-kB- 
regulated genes reported to be TNFa-responsive in mouse 3 T3 fibroblast 
cells (NFKBIA, CASP4, CCL5, TNFAIP3, CCL2, ZFP36, RIPK2, TNFSF10, NFKBIE, 
IL6, CCL20, ICAM1, TNFRSF1 A, TNFRSF1B, SQSTM1, NRG1, SOD1, IL1RL1, 
HIF1A, ERBB2) [19]. Error bars represent the 95 % confidence interval as 
determined by the Uncertainty statistic. Scores that failed the Specificity 
criterion (Specificity p-value > 0.05; 1000 comparable HYPs) are shaded in 
gray. 

Additional file 5: The aggregated IKK/NF-kB signaling HYP. Each 
row of the table contains a causal statement extracted from the Selventa 
Knowledgebase and obtained from the aggregation of the individual 
HYPs of the IKK/NF-kB signaling causal network model (Additional file 6). 

Additional file 6: The IKK/NF-kB signaling causal network model. 

The full causal model is given (top), along with a schematic of the basic 
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